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Abstract 

A two-parameter trial condensate wave function is used to find an approxi- 
mate variational solution to the Gross-Pitaevskii equation for Nq condensed 
bosons in an isotropic harmonic trap with oscillator length do and interacting 
through a repulsive two-body scattering length a > 0. The dimensionless pa- 
rameter A/q = Noa/do characterizes the effect of the interparticle interactions, 
with A/q ^ 1 for an ideal gas and Ao » 1 for a strongly interacting system 
(the Thomas- Fermi limit). The trial function interpolates smoothly between 
these two limits, and the three separate contributions (kinetic energy, trap po- 
tential energy, and two-body interaction energy) to the variational condensate 
energy and the condensate chemical potential are determined parametrically 
for any value of Ao, along with illustrative numerical values. The straightfor- 
ward generalization to an anisotropic harmonic trap is considered briefly. 
PACS numbers: 03.75.Fi, 05.30.Jp, 32.80.Pj, 67.40.Db 



I. INTRODUCTION 



The remarkable observation of Bose- Einstein condensation of dilute ultracold alkali atoms 
in a harmonic trap has stimulated much theoretical activity, based mostly on the 

Gross-Pitaevskii (GP) equation that is believed to describe accurately the experimental 
situation well below the onset of Bose-Einstein condensation. In this picture, the two-body 
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interaction is characterized by the s-wave scattering length a (typically, a ~ 10 nm and 
positive, but a is negative for "^Li [Q). The actual traps are anisotropic, but, for simplicity, the 
present work concentrates on an isotropic trap with frequency and oscillator length do = 
^h/ mujQ (typically, do ~ a few yum). Although the nonlinear GP equation is readily solved 
numerically for the selfconsistent wave function \1' of A^o condensed atoms in an isotropic 
trap with a total of atoms, it is often useful to have an analytic approximation 
that provides a semiquantitative description of the condensate density no(r) associated with 
the Nq < N condensed atoms. For example, Baym and Pethick ^ proposed two very 
different approximations: a Gaussian trial function in the limit of a small condensate (A/q = 
No\a\/do <^ 1) and a Thomas-Fermi (TF) approximation in the limit of a large condensate 
(A/q = NQ\a\/do ^ 1). The present work studies a simple variational trial function that 
includes the two limiting cases and yields a useful interpolation for intermediate values of 
Afo- 

The basic formalism is reviewed in Sec. II, and the explicit variational trial function 
is analyzed in Sec. Ill for positive scattering length a > 0. Typical numerical values are 
presented in Sec. IV, along with the various components of the energy per particle and the 
chemical potential. The quite different situation for negative scattering length is discussed 
briefly in Sec. V, and the generalization to an anisotropic harmonic trap is considered in 
Sec. VI. 



II. BASIC FORMALISM 

In the special case of a uniform dilute Bose gas at T = K subject to periodic boundary 
conditions, Bogoliubov [Q made the essential observation that the total number of condensed 
atoms is large {Nq ^ 1). More generally, for a nonuniform condensate, the occupation of 
the various single-particle states can be characterized by creation and annihilation operators 
a] and aj that obey Bose commutation relations [aj, al] = 6jk- The condition A'o ^ 1 implies 
that the operators for the condensate mode can be treated simply as numbers PJI(]|], because 
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their commutator [oq, Oq] = 1 is small compared to their value ~ aj) ~ \/Nq. As a result, 
the condensate is described by a conventional wave function \E'(r) and a "thermodynamic- 
potential" functional [obtained as a Legendre transformation from the "energy" functional 
EoiNo)] 

Koii^o) =Eo- fioNo = j Sr-^* {T + U - /io) ^ + 2rfm-^ j d\ |^|^ (1) 

that gives the dominant contribution for small noncondensate number N' = N — Nq <^ N. 
Here Eq(Nq) is the condensate energy, T = — /i^V^/2m is the kinetic energy and U{r) = 
^muj^r'^ is the isotropic harmonic trap potential. The condensate wave function is normalized 
to the total condensate number 

d^r\^\'^ = No, (2) 

so that the condensate particle density is simply no(r) = |^(r)p; the Euler-Lagrange equa- 
tion for the functional in Eq. (P is just the nonlinear Gross-Pitaevskii (GP) eigenvalue 
equation 

(T + U + V)^ = fio"^, (3) 

where V^(r) = ATcah'^no(r)/m is the mean (pseudo)potential energy per condensed particle. 
In the present case of a harmonic trap, /io must be adjusted to ensure that the wave function 
obeys the asymptotic boundary condition that \1' ^ as r ^ oo. In addition, /xq obeys the 
thermodynamic relation 

or, equivalent ly, 

N„ = (5) 

If the exact condensate wave function is known, then the exact condensate energy is given 
by 



Eo = {T) + {U) + l{V), (6) 

where the angular brackets denote an expectation value evaluated with the condensate wave 
function (namely, (■■■) = / d^r \E'* ■ ■ ■ \E'). A combination of Eqs. (^) and gives 

^,No = {T) + {U) + {V), (7) 

because V itself is proportional to Nq [this last result also follows directly from the GP 
Eq. m- 



III. VARIATIONAL TRIAL FUNCTION 

It is convenient to introduce dimensionless units, with do and huo as the length and 
energy scales. In this case, the dimensionless kinetic energy and trap potential energy 
become T = — and U = |r^. As a simple variational trial function for the isotropic 
condensate, take 

^(r) = Co (1 - rVi?')(^+^)/2 e{R - r) (8) 

that depends on the two dimensionless parameters A and R. The normalization constant Cq 
follows immediately from Eq. @ 

° 27ri?3r(|)r(2 + A)" 

For A — » 0, the trial function obviously reduces to the TP approximation \E'T'^(r) oc 
(1 — r"^ /R^)^^'^ 6(R — r), which follows by neglecting the kinetic energy relative to the trap 
potential energy and the interparticle potential energy; in this case U — fio + V vanishes 
in the region where \1/ 7^ 0, and /iq ~ |-R^ for an isotropic harmonic trap. The condensate 
density in this TP limit is parabolic, with no(r) oc (1 — r^/R"^) 6{R — r), showing that R is 
the (large) dimensionless TP condensate radius, and it is easy to verify that i?^ ~ 15A/o, 
where A/q = N^a/dQ ^ 1, and that Eq/Nq ^ |/io. 
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The opposite case occurs for A — oo. In this hmit, a detailed analysis (see below) shows 
that A/q ~ V327r/ A <^ 1, and that R ~ y/X ^ 1; in addition, a standard asymptotic identity 



involving the Gamma function [11 



^^^^z'^-'[l + Oiz-% for|.|->oo, (10) 

verifies that the trial function reduces to a normalized isotropic Gaussian with condensate 
density no{r) = |^g('^)P = Noir'^^"^ limA^oo (1 - r^/X)^^^ ~ iVovr"^/^ exp{-r^), as used 
by Baym and Pethick |^ (in this low-density, near-ideal, limit, the condensate density is 
independent of R, and the mean-square condensate radius is |). 

For general A and R, the various expectation values are readily evaluated: 

(T) = f AroA(A)i?-2, where A(A) = (1 + A)(f + A)/A; (11) 

(U) = I A^o B{X) i?^ where 5(A) = (| + A)-^ (12) 

[note that the mean-square radius (r^) is just 2{U)/Nq = |-B(A) R^, which depends on both 
A and R]; and 



r(3 + 2A) 



{V) = 2NoAfoC{X)R-', where C(A) = -^^^-^^^ ( . (13) 



r(|)r(| + 2A) 

Thus, the variational energy functional becomes 



r(| + A) 



r(2 + A) 



EpjR, A) _ 3 

No ' 4 



A{X)R-^ + ^,B{X)R^+AfoC{X)R-\ (14) 



Since A, B, and C are all positive, this expression obviously has a minimum for some positive 
R, and the derivative with respect to R immediately gives the condition 

- AR + B R^ = 2AfoC. (15) 

In addition, the derivative with respect to A yields the second condition 

I A'iX) R-^ + I B'iX) R^ + Xo C'iX) R~^ = 0, (16) 
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where C'{X)/C{X) is a linear combination of if) functions [0. Eliminating Mq/B? from 



Eqs. ([T5| ) and (|I6|) leads to an explicit equation for i? as a function of A: 

o4 _ ^ C'/C - I A' /A 

BC'/C+\B'/B' ^ ' 

and Eq. (^) then gives the corresponding value of Mq. 

A. Large condensate (A/q » 1 and A ^ 1) 

As shown in Sec. IV, these expressions are readily evaluated for any positive A, but it 
is helpful first to study two limiting cases, where analytic expressions can be obtained. For 
A ^ 1 and Aq ^ 1, a detailed expansion yields 

C'/C ^ Co + Ci A + ■ ■ ■ , (18) 

where Cq = f and Ci = -fyr^ + f + A + i| ~ -0.04432, so that C'/C + \B' / B ^ 
(Ci + ^) A + ■ ■ ■ . It is then easy to obtain the explicit result 

10^ 

^'- 8(C. + j|)A3 - " 3.600 A-/', (19) 

along with no{r) = |\1>tf('^)P ~ y (No/ttR^) (1 — r"^ / R^) and the specific results /io ~ | R^, 
Eq/Nq ^ n ^^'^ A/q ~ ^ -R^ given in the paragraph below Eq. (||). 

B. Small condensate (TVq < 1 and A > 1) 

In the opposite limit A 3> 1 and TVq ^ 1, it is simpler to return to Eq. (|15|), whose 

solution is of the form R^Mq) ^ Rq + {C/2A)N'q H with R^ = B/A. To first order in 

Ao ^ 1, the condensate energy is simply 

Eo(A)_i.o[i?o(A),A]_3 (20) 



Ao Ao 2 V ^ ^ ^ ^ i?o(A)3 

A direct expansion shows that A{X)B{X) ^ 1 + fA^^ H , i?o(A) ^ a/A (1 + |A"^ H ), 

and C(A) ^ A3/2 (27r)-V2 (i + 69^-1 + ...). thus, 
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+ 



No 2 8A2 V 16 A 

This expression clearly has a minimum as a function of A, occurring at 



Mo 



IT' 



the corresponding condensate energy 



(21) 



(22) 



(23) 



precisely reproduces that found directly for the isotropic normalized gaussian trial function 



IV. NUMERICAL RESULTS 

For any positive value of the variational parameter A, it is straightforward to evaluate the 
other dimensionless variational parameter R from Eq. {^i7\), the interaction parameter A/o = 
Noa/do that specifies the condensate number from Eq. (^5[), and the various contributions 
to the condensate energy per particle and the condensate chemical potential. Typical results 
are shown in Table 1, illustrating the continuous interpolation between the TP limit (large 
condensate with A/q ^ 1 and A <^ 1) and the Gaussian limit (small condensate with A/q ^ 1 
and A ^ 1). In the TP limit, the small-A approximation in Eq. ([T9|) yields R ~ 113.8 for 
A = 0.01, with the following TP values A^o ^ j^R^ ~ 1.275 x 10^ and ^ ^ 6480; 
these values are close to the respective variational ones 114.4, 1.297x10^, and 6526 given in 
Table 1. Similarly, the Gaussian (large-A) limit takes R ~ \/A; for A = 100, this relation 
gives R = 10, and the approximate expressions in Eqs. (|2^ ) and (p3D yield Ao ~ 0.100 and 
Eq/Nq ^ 1.540, which should be compared to the respective variational values 0.105 and 
1.541 from Table 1 (which also yields the mean-square radius ~ 1.56, close to that for the 
ideal gas). Note that the physical parameter Ao = Noa/do is a single- valued function of A 
that decreases monotonically with increasing A. 
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As an example of the utility of this variational approach, consider the situation for 
A = 1, where the ratio {T)/{Eo) ^ 0.022 is already small (so that /iq ^ (f/ + V)/No), yet 
the variational condensate density oc (1 — r^/i?^)^ differs greatly from the parabolic 

TF condensate density |\1/tfP oc 1 — r'^/R^. Thus the elementary criterion {T)/{Eq) ^ 1 
alone is insufficient to ensure that the condensate density approaches the TF form. 

It is important to recall that the variational method provides only an upper bound for 
the total condensate energy, so that the separate contributions are not necessarily accurate. 
In particular, Dalfovo, Pitaevskii, and Stringari have used a boundary-layer formalism 



to show that the kinetic energy per particle in the TF limit is oc i? ^ In i?, which differs 
significantly from that found variationally. 

V. BEHAVIOR FOR NEGATIVE SCATTERING LENGTH 

The stability of a Bose condensate with negative scattering length has become of interest 
in connection with the experimental study of ^Li . It is well known that a uniform Bose 
condensate is unstable for arbitrarily small negative interactions (the speed of sound be- 
comes imaginary ||9|Jl0|l); this behavior is somewhat reminiscent of Jeans's long- wavelength 
(gravitational) instability for an infinite uniform stationary medium fl^, where the disper- 



sion relation is simply that for a bulk plasma oscillation with the sign of the interaction 
reversed. The external confining trap completely alters the situation, however, for it elim- 
inates wavelengths longer than the characteristic dimension of the trap (like electromag- 
netic standing waves in a cavity). Indeed, the Bogoliubov dispersion relation for a uniform 
medium with negative scattering length —\a\ becomes |]9|jl0[| El = (h^k/2m)'^ (A;^ — 167r|a|no), 



which is stable only for wavenumbers greater than kc = yl67r|a|no. For an order-of- 
magnitude estimate, take uniform density and /cmin ~ TT/2do (appropriate for a spherical 
square well of radius do); the condition fcmin > yields the approximate stability criterion 
\J\fo\ < A/qc = NQ\a\/do = 7r^/48 ~ 0.206. For comparison, a numerical study of the GP 
equation for negative scattering length in an isotropic harmonic trap M yields the exact 
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critical value A/oc = 0.573. 

In the present context, it is interesting to examine how the variational calculation treats 
a negative scattering length. Since the trap can remain stable at T = K only for relatively 
small condensates, it is appropriate to consider the behavior for small negative values of 
Aq. Equation (^Tj) remains valid, but the altered sign of Aq means that Eq{X) achieves 
its minimum only for A — oo, corresponding to a Gaussian trial function. This particular 
situation has already been examined in |T^Jl5[l , where the local minimum in the variational 
energy has been shown to disappear at the critical value |A/oc| = (87r/25v^)"'^/^ ~ 0.671, 
somewhat larger than the exact value found numerically in 0. 



VI. ANISOTROPIC HARMONIC TRAP 

The variational approach has the appealing feature that it is readily generalized to an 
anisotropic harmonic trap potential, with 

f/(r) = lmE, ^|r|, (24) 

where j = 1, 2, and 3 runs over the three coordinate axes. Correspondingly, the three (in 
general, distinct) oscillator lengths are given by (i| = h/muj. As a natural trial function, 
take 

*w -0 (i - 1: -g) Ki-Esm)- p«) 

where {Rj} is a set of three dimensionless variational parameters. The normalization con- 
stant is easily evaluated with Eq. (||) by going to "scaled" variables Xj = (rj/djEj), yielding 
the generalization of Eq. (|^) 

= ^0^(1 + ^) .26) 

2nR,R,RsTil)T{2 + Xy ^^^^ 

The three relevant expectation values are also readily found for this more general trial 
function 
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{T) = \noA{X)J:^, (27) 



{U) = ^NoB{X)J2hu;,Rl (28) 



and 



{V)=2N,nu,§£^, (29) 

where ciJ^ = iUiuj2^3, A/q = N^a/ {did2d^Y^^ ^ and the functions A, 5, and C are the same as 
for an isotropic trap, given in Eqs. (|ll])-(|l3D. Thus, the general variational energy functional 



depends on the four variational parameters -Ri, i?2, -Rs, and A. 

It is not difficult to minimize this functional with respect to each of the first three 
parameters. The three conditions dEo/dRj = yield the equations 

- A/R^j + BR] = 2Duj/ujj, (31) 

where D = JVq C {X) / R1R2R3 is the same in each equation. This set of conditions can be 
solved formally to express Rj as 

uJD 



R 



'•5 , 1 



+ \-ls]^ (32) 



j,B \ B \u,B 

where the sign of the root is chosen to give the correct limit for a ^ (the ideal gas). Since 
-D^ = {JVoCy / RIR2R3, these expressions immediately provide one relation between D and 
A. In addition, the minimum condition dEo/dX = can be combined with Eqs. (^) and (|32D 
to provide a second relation between D and A, thus determining both quantities, as well as 
the three dimensionless parameters Rj and Aq. Edwards et al. |[T^ and Dalfovo and Stringari 



have used distinct numerical procedures to evaluate the ground-state condensate wave 
function for anisotropic but axisymmetric traps, using the physical parameters appropriate 
to [|T| ; it would be interesting to compare the present variational approach with their explicit 
results, but such a numerical study remains to be performed. 
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TABLES 

TABLE I. Typical variational physical parameters that characterize a spherical condensate in 
an isotropic harmonic trap 



A 


R 


Mo 




(C/)/Ao = i(r2) 




Eo/No 




0.01 


114.407 


1.2974x10*^ 


0.0145 


2796.77 


1864.51 


4661.29 


6525.80 


0.1 


21.2428 


2.685x10^ 


0.0475 


94.0115 


62.6426 


156.702 


219.3 44 


1 


5.3967 


153.92 


0.1803 


4.8541 


3.1159 


8.1503 


11.2662 


10 


4.3724 


1.5756 


0.5394 


1.0621 


0.3484 


1.9499 


2.2984 


100 


10.3758 


0.1054 


0.7212 


0.7801 


0.0393 


1.5406 


1.5799 
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